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^ : ABSTRACT 

^: 

. We model how repeated supernova explosions in higli-redsliift dwarf star- 

>• . burst galaxies drive superbubbles and winds out of the galaxies. We compute 

^ . the efficiencies of metal and mass ejection and energy transport from the galactic 

. potentials, including the effect of cosmological infall of external gas. The star- 

burst bubbles quickly blow out of small, high-redshift, galactic disks, but must 
compete with the ram pressure of the infalling gas to escape into intergalactic 
space. We show that the assumed efficiency of the star formation rate domi- 
nates the bubble evolution and the metal, mass, and energy feedback efficiencies. 
With star formation efficiency = 0.01, the ram pressure of infall can confine 
the bubbles around high-redshift dwarf galaxies with circular velocities Vc ^ 52 
km s~^. We can expect high metal and mass ejection efficiencies, and moderate 
energy transport efficiencies in halos with fc ~ 30 — 50 km s~^ and /=„ ~ 0.01 
as well as in halos with Vc ~ 100 km s~^ and /* ^ 0.01. Such haloes collapse 
successively from l-2a peaks in ACDM Gaussian density perturbations as time 
progresses. These dwarf galaxies can probably enrich low and high-density re- 
gions of intergalactic space with metals to 10~'^-10~^ Zq as they collapse at 2; ~ 8 
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and z < 5 respectively. They also may be able to provide adequate turbulent 
energy to prevent the collapse of other nearby halos, as well as to significantly 
broaden Lyman-a absorption lines to v„ns ~ 20-40 km s~^. We compute the 
timescales for the next starbursts if gas freely falls back after a starburst, and 
find that, for star formation efficiencies as low as < 0.01, the next starburst 
should occur in less than half the Hubble time at the collapse redshift. This 
suggests that episodic star formation may be ubiquitous in dwarf galaxies. 

Subject headings: hydrodynamics — shock waves — galaxies: dwarfs — galaxies: 
starburst — galaxies: high-redshift — galaxies: evolution 

1. Introduction 

1.1. Overcooling and Angular Momentum Problems 

The theory of hierarchical structure formation in a universe dominated by cold dark 
matter (CDM) predicts that galaxies assembled from mergers of dark matter halos that 
gravitationally collapsed from primordial ffuctuations. Within the same framework, the first 
luminous structures to form have sub-galactic masses, given by the Jeans mass of the gas. 
This in turn is determined by the cooling mechanisms dominant during each cosmological 
era, since CDM models do not define a lower limit to the scale of inhomogeneities in the 
overall mass distribution. Protogalactic gas clouds of mass Mh ~ 10^ M© begin to collapse 
as ~ 3(7 peaks around redshift z ~ 10 and become dwarf galaxies. 

The hierarchical scenario predicts the formation of numerous dwarf-sized halos at high 
redshift. In these halos, the dissipative collapse of gas is very efficient, since the cooling time 
T oc p"^ oc {l + z)~^, where p is the mass density. This is often referred to as the "overcooling 
problem" (White & Rees 1978; White & Prenk 1991). Overcooling leads to a prediction of 
far more dwarf galaxies than observed, as well as to an angular momentum problem. When 
dense, cooled, protogalactic gas clouds or dwarf galaxies interact within a larger system such 
as the halo of a present-day disk galaxy, dynamical friction transfers their orbital angular 
momentum to the surrounding dark matter halo. The result is a galaxy with too much mass, 
rotating too slowly, with angular momentum deficient by a factor of ~ 25 relative to that 
observed (e.g. Katz & Gunn 1991; Navarro & Benz 1991; Navarro & White 1994). 

Both of these problems can be resolved if a mechanism can be found to keep the gas 
diffuse until the peak of present-day disk (L^,) galaxy formation. Supernova feedback from 
young stars is one commonly invoked mechanism to do this. A number of cosmological 
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hydrodynamic simulations have demonstrated that stellar feedback plays a major role in the 
formation of galaxies (Katz 1992; Navarro & White 1993; Yepes et al. 1997; Gerritsen & 
Icke 1997; Hultman & Pharasyn 1999; Scannapieco et al. 2001; Thacker & Couchman 2001; 
Sommer-Larsen et al. 2002; Springel & Hernquist 2003). However, the qualitative outcomes 
of simulations vary drastically with the feedback scheme applied. Recently, smoothed particle 
hydrodynamic (SPH) simulations by Scannapieco et al. (2001) and Springel & Hernquist 
(2003) show that overcooling in low- mass halos is significantly suppressed by the galactic 
outflows that they model, since the outflows slow down the accretion of gas in the halos and 
also strip baryons from neighboring halos ("baryonic stripping"). The latest high- resolution 
(~ 10^/i~^MQ per gas particle) SPH simulations by Thacker and Couchman (2001) and 
Sommer-Larsen et al. (2002) show that their treatment of thermal feedback with suppressed 
cooling succeeds in reproducing a present-day disk galaxy with angular momentum that is 
deficient only by a factor of a few. 

However, these cosmological simulations still fail by orders of magnitude to resolve the 
physics of star formation and feedback. The coupling of feedback energy with interstellar gas 
through shocks, the formation of galactic outflows, and the coupling of the outflow energy 
with halo gas are not physically represented, especially in SPH simulations with resolution 
at the galactic scale far too low to accurately resolve shocks. Therefore, it is important to 
model the collective action of multiple supernovae and study the formation of superbubbles 
and galactic winds in single dwarf galaxies, using high-resolution, hydrodynamic simulations, 
in order to compute the feedback efficiencies (e.g. Mac Low & Ferrara 1999; hereafter MF99). 
Such models must also include cosmological infall of gas at high redshift, because the ram 
pressure of the infall may influence the evolution of the superbubbles and the galactic winds. 



1.2. Intergalactic Metals 

Observations of Lya-absorbing clouds reveal the presence of metals. These clouds are 
thought to be regions of enhanced intergalactic medium (IGM) density distant from proto- 
galaxies (see review by Rauch 1998). To transport metals into these regions, some mechanism 
such as supernova- driven galactic outflows must act. The observations show that Lyman for- 
est clouds with neutral hydrogen column density with log N{H I) > 14 are metal enriched 
to Z ^ IQ-'^ - IQ-^ Zq (Cowie et al. 1995; SongaUa & Cowie 1996; Eflison et al. 2000), and 
that the mctallicity remains roughly constant throughout the redshift range 1.5 < z < 5.5 
(Songaila 2001). It is also interesting to note that log N{Hr) ^ lA at z ^ 3 marks the 
transition between continuous fllamentary structures and voids (e.g. Zhang et al. 1998). 

The main candidates for the polluters of the IGM are starburst dwarf galaxies. The 
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absence of turbulent motions observed in low-density Lyman a clouds at 2; ~ 3 (Ranch et 
al. 2001b) suggests that metal enrichment was completed at very early times (2; ^ 5), when 
the physical volume of the universe was smaller, and contained the numerous dwarf galaxies 
predicted by the CDM model of galaxy formation (Ranch et al. 2001b; Madau, Ferrara, & 
Rees 2001; Scannapieco, Ferrara, & Madau 2002). On the other hand, the median Doppler 
parameters are significantly larger in the Lyman a forest than predicted by cosmological 
simulations, suggesting additional energy injection provided by late Hell reionization or 
supernova-driven winds (Meiksin, Bryan, & Machacek 2001). In CiV systems with higher 
column densities, a substantial velocity scatter over length scales of a few hundred parsecs is 
observed, however, suggesting they have been influenced by galactic feedback more recently 
(Ranch et al. 2001a; see Ranch 2002 for summary). Galactic outflows both at very high 
redshift and at lower redshift 2; ~ 3 seem to play signiflcant roles in enriching various regions 
of the IGM. 

Cosmological sinmlations have suggested that the IGM at ^ ~ 3 can be metal-enriched 
to Z = 10^^-10^^ Zq by merging of or by outflows from dwarf galaxies (Gnedin 1998; 
Theuns et al. 2002; Thacker et al. 2002) without dynamically disturbing the observed, low- 
density, Lya clouds. These simulations have the same problem as earlier that the feedback 
schemes differ based on different assumptions about the coupling of supernova energy to the 
surroundings, and the simulations fail badly to resolve star formation and the dynamics of 
bubbles. 



1.3. Starburst Wind Models 

The effects of repeated supernova explosions from starbursts appear to be a central piece 
of physics necessary for understanding the role of stellar feedback in both galaxy formation 
and IGM pollution. The effects of starbursts from dwarf galaxies on the surrounding inter- 
stellar medium have been studied in the past both analytically and numerically (Mathews & 
Baker 1971; Larson 1974; Saito 1979; Dekel & Silk 1986; De Young & Heckman 1994; Silich 
& Tenorio-Tagle 1998; MF99; D'Ercole & Brighenti 1999). For example, MF99 studied the 
effects of repeated supernovae on local dwarf galaxies with gas mass Mg = 10^-10^ M0 using 
hydrodynamic simulations, and made a parameter study of mass and metal ejection efficien- 
cies. The main results of MF99 are that mass loss is very inefficient except in the lowest 
mass halos, while a substantial fraction of the hot, metal-enriched gas escapes from any of 
the potentials they studied. These results were confirmed by D'Ercole & Brighenti (1999) 
with a better treatment of thermal conduction. Recently, Mori, Ferrara, & Madau (2002) 
studied the same problem, but in a spherical, non-rotating, galaxy embedded in a halo with 
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Mh — 10^h~^ Mq at 2; = 9. They find that, with star formation efficiency /* = 0.1, about 
30% of the available supernova energy is transferred to the surrounding gas as kinetic energy. 
The starburst bubbles expand over ^ lOR^, where is the virial radius of the halo, com- 
parable to the mean proper distance between neighboring low-mass systems. Their results 
support the suggestion that IGM metal enrichment occurs early. The recent study by Wada 
& Venkatesan (2003) improves on Mori et al. (2002) by precalculating a thin, dense disk with 
an inhomogeneous ISM in a given halo potential and by including self-gravity. They perform 
two three-dimensional models of a 10^ Mq galaxy with star formation efficiency /* = 0.14 
and 0.014 in our terms. 

We follow the work of MF99 and extend their study to high redshift. Although the 
approachs taken by Mori, Ferrara, & Madau (2002) and Wada & Venkatesan (2003) are more 
realistic in some ways, the former neglects rotational flattening of the galaxy, and both are 
far more expensive computationally. The MF99 approach enables us to do a parameter study 
on the feedback process by exploring a wide range of galaxy masses, formation redshifts, and 
star formation efficiencies. However, note that our study is limited to axisymmetric geometry 
and also to single star formation sites at the centers of disks with a smooth ISM. 

We perform careful, high-resolution models of the result of multiple supernova explo- 
sions in single dwarf galaxies at high redshift, using a hydrodynamic code, ZEUS-3D (Stone 
& Norman 1992; Clarke 1994), and compute the efficiency of metal and mass ejections and 
energy transport from the galactic potentials. We include the evolving dark matter halo 
potentials, and the cosmological infall of halo gas computed with a one-dimensional hydro- 
dynamic code (Meiksin 1994). The bubbles in MF99's local galaxies were solely governed by 
their interaction with the interstellar medium (ISM), and so high metal ejection efficiencies 
were guaranteed once the bubbles blew out of the galactic disks. However, the bubbles in 
our high-redshift galaxies must still fight the ram pressure of the infalling halo gas after 
blowout before the metals, swept-up mass, and energy escape to intergalactic space. We also 
compute the accretion timescales for any mass swept up by the bubbles, but bound by the 
potential to fall back to the center. 

We model systems that we think hkely to host a starburst. It is not certain whether 
first generation systems were efficient in ejecting metals and transporting energy to the IGM, 
because an early cosmic UV background suppreses the formation of stars inside the systems 
by photodissociating their only cooling agent, molecular hydrogen (Haiman, Abel, & Rees 
2000), or by raising the entropy floor (Oh & Haiman 2003) However, later systems that 
collapsed due to atomic hydrogen cooling may have been more robust allowing starbursts to 
occur (Madau, Ferrara, & Rees 2001; Scannapieco, Ferrara, & Madau 2002; Oh & Haiman 
2002). 
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We choose seven of these second generation systems, with 5 x 10^ Mq < Mh < 5 x 10^° 
Mq at 3 < z < 13 with star formation efficiencies /* = 0.001, 0.01, and 0.1. Our goal is 
to develop a generally useful description of supernova feedback in terms of metal ejection 
efficiencies C,metai, mass ejection efficiencies ^, and energy transport efficiencies that can 
be widely used in large-scale cosmological simulations for the study of galaxy formation and 
metal enrichment. However, it will be shown in this paper that parameterization of feedback 
is very difficult, and that we eventually need a higher-resolution, adaptive mesh refinement, 
cosmological simulation with a realistic treatment of star formation to do so as reliably as 
we might wish. Nonetheless, we are able to reach several key conclusions with the present 
set of computations. 

The paper is organized as follows. In §2 we analytically predict the effects of ram pressure 
on the confinement of the starburst bubbles. In §3 we discuss our numerical methods, 
including our models of cosmological infall (§3.1), starburst blowout (§3.2), and the tracer 
field we use to follow metal-enriched gas in the blowout model (§3.3). In §4 we describe 
our choices for galaxy size (§4.1), disk model (§4.2), and star formation feedback (§4.3). We 
give our results in §5: metal and mass ejection, and energy transport efficiencies in §5.1, the 
effects of external pressure on the bubble evolution in §5.2, and the timescales for material 
to fall back in §5.3. Our conclusions follow in §6. Throughout the paper, we apply ACDM 
cosmology with Qq — 0.37, fl\ = 0.63, fib = 0.05, h — 0.7, as ~ 0.8 and a primordial 
composition of H:He=12:l. 



2. Pressure Confinement 

We first try to estimate analytically how the ram pressure of infalling gas affects the 
evolution of a starburst bubble, in order to estimate the ranges of halo mass Mh, redshift 
z, and star formation efficiency /* over which the ram pressure confinement of a bubble 
suppresses feedback. 

In this section, we consider a bubble in a power-law halo with gas density p = po^ 
Isothermal halos with constant angular velocity have a = 2. We take the bubble to be 
driven by a mechanical luminosity Lm for the 50 Myr lifetime of the least massive star able 
to go supernova after an instantaneous starburst, and then to evolve freely thereafter. A 
bubble big enough to interact with the halo will blow out of the highly stratified, exponential 
atmosphere of a high-redshift galactic disk very quickly. After it blows out, its evolution is 
largely determined by the infalling halo gas. Since the bubble accelerates as it blows out 
of the exponential ISM (Kompaneets 1960), it is difficult to exactly predict its expansion. 
Therefore, we make the simplifying assumption that the bubble has reached the power-law 
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halo by 50 Myr, which is supported by our detailed numerical calculations described in 
section 5 below. Thereafter, we compare the interior pressure of the bubble with the ram 
pressure of the infalling halo gas. 

After the starburst ends at a time tf — 50 Myr, the bubble continues to expand due to 
the residual energy of supernovae stored in the hot interior gas and in the shell. However, 
we do not know how much of the injected energy L^tf remains in the interior, rather than 
being dissipated by radiative cooling enhanced by shell fragments mixing with the interior 
gas, and how much of it is transferred to the kinetic and thermal energy of the swept-up 
shell. Therefore, we assume a fraction u of the injected energy Lmtf drives the bubble after 

tf- 

From the Rankine-Hugoniot jump conditions, we have 

P - (1) 

where P is the pressure behind the shock front, p — PqR~^, the radius of the shock is 
and Ui is the velocity ahead of the shock, in a frame of reference traveling with the shock 
velocity Rg. The preshock velocity is 

Ui^ Rs + Vin{Rs), (2) 

where the infall velocity of gas is Vin{Rs)- We assume that the pressure P behind the shock 
equals the mean pressure P of the gas within the spherical volume enclosed by Rg, and a 
fraction z/ of the total injected energy L^tf gives the pressure (Spitzer 1978). (It is equivalent 
to setting = 1 in the Kompancets approximation; however, in reality, it is not always true 
that P = P: for example, P ~ 2P in a Sedov blast wave solution [Sedov 1959]). We can 
then derive the shock velocity, 

Rs = 0.46 (^^) Ri^-^y' - vURs), (3) 

where we have taken 7 = 5/3 for an adiabatic bubble shock (but the coefficient varies by only 
15% for a radiative shock with 7 = 1), and we assume the thermal pressure of the IGM to be 
much smaller than its ram pressure. The shock velocity driven by the mechanical luminosity 
prior to tf exceeds that given by equation (3), but the bubble quickly decelerates without 
pressure support from the central energy source. 

We approximate the infall velocity of gas. 



Vin{r) - {2GMh/r) 



(4) 
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as that of a free fall onto a point mass with the mass Mh of the halo, as long as the radius 
r exceeds the virial radius Ry (given below by Eq. [16]). We completely neglect the angular 
momentum of the infalling gas, so Vj„ is an upper limit. The bubble will stall when Rg — 
by ram pressure exerted by the infalling gas. Then, we can derive 

^,„(i?,)>o.46(^^y^'i?^3)/^ (5) 

as the criterion for ram-pressure confinement. Both sides of the equation are functions of 
the bubble radius Rs, but we can eliminate the radial dependence in an isothermal halo 
with a — 2 and compare both sides of the equation directly without knowing where the 
bubble is. Note that the left side of the equation is just under the escape velocity of a halo 
Vesc ~ {GMh/ R^y/'^ a,t r > Rv, and the right side is the shock velocity of a bubble if there 
is no infalling gas. 

The coefficient po for the gas density of a spherically symmetric, power-law medium is 

when a given total halo mass Mh is found within the virial radius R^. The mechanical 
luminosity L„ can be expressed in terms of star formation efficiency as 

- (3.6 X 10- e.g (^) (4) , (T) 

based on the Starburst 99 model (Leitherer et al. 1999) described in §4. 

Substituting equations (4), (6), and (7) into equation (5) with a = 2, we find that 
external ram pressure can confine a bubble at a given redshift z when the total halo mass 

(n \ 3/2 / 1 I \ -3/2 
d^) (^) ■ («' 

We plot the above equation, the minimum halo mass for pressure confinement for different 
values of i//* in Figure 1. The starburst bubbles in the dwarf halos that we select for our 
study (filled diamonds: see §4.1) are predicted to freely escape to intergalactic space if 
/* > 0.1 and to be completely confined by ram pressure of the infalling gas if /* < 0.001 and 
v < 0.5. We discuss Figure 1 further in §4.1. For a given i/f^, both the minimum halo mass 
for pressure confinement in equation (8) and the halo mass with a fixed virial temperature 
Ty or circular velocity Vc have the same redshift dependence, Mh oc (1 + z)^^^"^. Therefore, 
we can select a group of halos in which the ram pressure of the infall can confine the bubbles 
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and so suppress the feedback efficiencies by a virial temperature or circular velocity cut-off: 
T* or V*. We can not expect high feedback efficiencies in dwarf galaxies with 

> T; ~ (1.0 X 10^ K) ' 



V, 



1 X 10-3 

> < ^ (52 km s-^) ( V^' . (9) 



1 X 10-3 



Without ram pressure from cosmological infall, a bubble will stall when its expansion 
is balanced by the thermal pressure of the IGM. If we follow the same procedure as above, 
we find that the bubble will only stall when its radius 

«.^(io«.)(r^)(^)". (10) 

where we have assumed an ionized IGM. 

The presence of ram pressure does generally appear to be important to the evolution of 
bubbles, and hence must be included when we compute the metal, mass, and energy feedback 
efficiencies in numerical simulations. 



3. Numerical Methods 

In this section, we outline the numerical methods we use to model feedback from dwarf 
galaxies. We first describe the model that we use for the evolution of the background halo, 
then the implementation of our models of starbursts. Finally we examine an important 
detail of our model: the performance of the tracer field that we use to trace hot, metal- rich 
material from bubble interiors. 



3.1. Spherical Halo Model 

We model evolving dark matter halo potentials and the cosmological infall of gas onto 
them using a method based on the one-dimensional hydrodynamic computations of Meiksin 
(1994), which were developed to study the structure and evolution of Lyman a clouds grav- 
itationally confined by dark matter minihalos. The infall code described by Meiksin (1994) 
included standard CDM cosmology; we modified it for a fiat, ACDM cosmology. We neglect 
the effects of the repulsive force exerted by dark energy as it is negligible at redshift z '>2. 
This can be seen by considering the expansion parameter 

H{tf = {a(t)/a{t)f = SttCpm/S + A/3 (11) 
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where A is the cosmological constant and pm — (1-88 x 10 "^^h^ g cm ^)Qo(l + As z 
increases, matter dominates. 

The infall code follows the evolution of spherical shells of dark matter based on a model 
of the spherical collapse of a Gaussian density perturbation (Bond et al. 1988). However, 
the dark matter shells are not allowed to interpenetrate as they collapse. Instead each shell 
is brought to rest at its virial radius, building a halo that reflects the core structure of the 
initial perturbation, so the dark matter halo profile will not perfectly reproduce the results 
of A^— body simulations, especially in the central regions. However, the central regions of the 
halo arc less important anyway, because the disk of the galaxy dominates the potential in 
that region. As we discuss in §4, we therefore include a separate disk term in the potential. 
The halo potential is more important for determining the density profile of the infalling gas 
at its edges, as well as capturing the details of the late evolution of the bubble as it chmbs 
out of the potential well. 

In the code, a hnear Gaussian perturbation is evolved starting at 2; = 200. For models 
that collapse at 2; > 5, we assume that the IGM is fully neutral with initial temperature 

= 550 K determined by the cosmic background radiation. For models that collapse 
at 2 < 5, we assume that H is ionized and He is singly ionized, with initial temperature 

Ti,g = 10^ K, and we set the minimum temperature to which the gas is allowed to cool to 
be Train = 10"^ K. This very roughly captures the effects of UV background radiation after 
reionization, which actually heats the IGM to ~ 2 x 10^ K (Sargent et al. 1980). 

We assumed that the reionization was completed at Zion ~ 6, and that the IGM was 
mostly neutral before Zion, based on the onset of Gunn-Peterson trough against background 
quasars (Djorgovski et al. 2001; Becker et al. 2001). However, the recent measurement of a 
large electron scattering optical depth by WMAP implies that a significant fraction of the 
IGM was already ionized by 2 ~ 17±5 (Kogut et al. 2003; Spcrgcl et al. 2003). If our models 
at 2 > 5 are evolved in the reheated IGM instead, most of the gas ends up being pressure 
supported in halos with 5 x 10^ K and Vc ~ 30 km s~^. However, the differences in 

density and pressure distributions appear only in central regions of the halos, where the gas 
cools to form disks. We discuss the effects of suppressed coohng in §4.2. 

The fluid equations are solved for baryons that cool by optically thin radiative recombi- 
nation losses from H and He, dielectric recombination losses from He, bremsstrahlung, and 
Compton cooling from scattering off the microwave background. Cooling is turned off in the 
central region where the potential is anyway unphysical to avoid thermal instability. CoUi- 
sional processes were included by Meiksin (1994) but were found to be unimportant outside 
of the accretion shock, so they have been neglected in our model to increase computational 
speed. We also neglect molecular hydrogen cooling. 
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The infall code produces a gas density distribution p oc r~° with a 2-2.3 for r < Ry 
and a ^2 outside, as predicted by the secondary infall model by Bertschingcr (1985). The 
density distribution at the center is not correctly modeled by the infall code. However, we 
replace the central gas distribution with a rotationally supported disk, as we describe in the 
next section. The lack of cooling in the central region does produce shocks that convert 
some infall energy to thermal energy, and so the distribution of thermal or ram pressure 
alone is not well represented physically in our infall solutions, but the total pressure remains 
accurate. The details of the energy and velocity profiles ultimately depend on how gas is 
removed by instabilities that make molecular clouds and stars, which we are still far from 
modeling. 

Although we neglect the large-scale filamentary web predicted by three-dimensional 
hydrodynamic simulations, the spherical infall code is probably not too bad for the scale of 
our interest, because virialized systems form at the intersection of the filaments, where the 
accretion is more spherical within a scale of order the virial radius closely corresponding 
to minihalos (e.g. Miralda-Escudc et al. 1996; Zhang et al. 1998). Adaptive mesh refinement 
simulations of first star formation by Abel et al. (1998; 2000; 2002) do show that gas accretion 
occurs along filaments, producing a highly structured, aspherical accretion shock even close 
to Ry. The results of our models in an idealized, spherical background can be interpreted as 
a lower limit to the effectiveness of metal, mass, and energy feedback, however, since a more 
dynamical, filamentary shape will always make it easier for galactic winds to find escape 
routes from dense regions. 



3.2. Starburst Models 

Our models of starburst blowout generally follow MF99 (see their § 4) and Fujita et al. 
(2003: see their § 4), so we briefly summarize our methods below. 

We compute the evolution of the blast wave from starburst supcrnovae with ZEUS, an 
Eulerian, finite-difference, astrophysical gas dynamics code (Stone & Norman 1992; Clarke 
1994), that uses second-order van Leer (1977) advection. We use a quadratic artificial vis- 
cosity to resolve shock fronts, and a linear artificial viscosity to prevent instabilities due 
to strong shocks in standing fiows. The former is necessary in any method that does not 
solve the Riemann problem, while the latter only acts at the zone-to-zone scale, and so does 
not affect larger-scale features of the solutions. We use the loop-level parallelized version 
ZEUS-3D, in its two-dimensional form. Runs were done on Silicon Graphics Origin 2000 
machines using eight processors, and typically took several days to a week, depending on 
the parameters of the galaxy. 



- 12 - 



We assume azimuthal symmetry around the rotational axis of the galaxy, and use a 
ratioed grid in both the radial and vertical directions to capture the multiple length scales of 
the problem. Our central, high-resolution grid extends to ~ 0.5i?„, and has typical zone sizes 
of a few parsecs at high redshift to a few tens of parsecs at lower redshift. Zone size increases 
by a factor of ten across our ratioed grid but in the outer parts of the grid still remains under 
a few tens to hundreds of parsecs respectively The grids always capture several Ry. The 
actual values for the different models are hsted in Table 1. The central resolution is high 
enough to resolve the stratified atmosphere of disks with more than 50 zones. This is found to 
yield converged results on the acceleration of bubbles, and is less than or at least equal to the 
typical size of a giant molecular cloud (~ 10-40 pc). We use reflecting boundary conditions 
along the symmetry axis and along the galaxy midplanc and infall boundary conditions on 
the other two axes. The infall boundaries are fed by interpolation from intermediate results 
of the spherical collapse model described in the previous section. 

To drive a constant luminosity wind, we add mass and energy to a source region with a 
radius of five zones and a volume V. We choose the input rates based on the analytic solution 
by Weaver et al. (1977) for a constant density medium with the density po and energy eg of 
the center of the galactic disk. The energy density in each cell of the source zone is increased 
at a rate e = Ljn/V, where is the mechanical luminosity from the starburst. The mass 
density is then increased at a rate p = poe/rieo to roughly keep the specific energy of the 
source region constant. 

The amount of mass in this wind accounts not only for supernova cjccta, but also mass 
that would have evaporated off the shells if we had included thermal conduction in our model. 
The evaporated mass does, in fact, dominate the mass within a realistic bubble (Weaver et 
al. 1977). Therefore, we add a larger amount of mass by choosing ri{L^) for a given L^, to 
reproduce the mass evaporation rate predicted by Weaver et al. (1977) for a spherical bubble. 
In this way, we keep the density and temperature within the superbubble fairly close to the 
analytic solution. See Fujita et al. (2003) for more details. With our choice of ri{Lm), the 
resulting interior temperature is always within a factor of two of the expected temperature. 

We use a cooling curve for primordial gas in coUisional ionization equilibrium (Suther- 
land & Dopita 1993) and also implement inverse Compton cooling appropriate for the high- 
redshift universe. We employ a sharp cutoff in the cooling at T = 10^ K. This not only 
distinguishes ionized gas which cools, from neutral gas which does not cool at high redshift, 
but also mimics the effect of UV background radiation at intermediate redshift, in the ab- 
sence of explicit models for photoionization and coUisional ionization in our code. In the 
high-redshift, dense universe, we expect bubble shell shocks to be nearly isothermal due to 
efficient coUisional cooling {Aune oc n^), so the coUisional ionization equilibrium cooling curve 
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is a reasonable approximation. Strictly speaking, the assumption breaks down at the shock 
fronts, where underionization will increase the cooling rate, and in the evaporative flow off 
the shell, where overionization will decrease the cooling rate. Using the equilibrium curve 
overestimates cooling in the interiors, ensuring that our feedback efficiencies are lower limits 
to the true value. 

The only place in the evolution of a bubble where substantial cooling occurs is in the 
dense, swept-up shell, where collisionally-excitcd radiative cooling determines the density. 
The cooling in the shell does not directly affect the dynamical evolution of a bubble, which 
is controlled by the hot interior. The cooling timescalc of the interior greatly exceeds the 
dynamical timescale of any realistic starburst bubble (Mac Low & McCray 1988). Instead, 
conductive and mechanical energy fluxes dominate its energy budget. Rayleigh- Taylor insta- 
bility occurs when the shell begins to accelerate during blowout, so that the dense swept-up 
shell of ISM and halo gas is supported by the low-density, hot gas of thermalized supernova 
ejecta and evaporated shell gas inside. 

The cooling curve is implemented in the energy equation in ZEUS-3D with a semi- 
implicit method. To solve the implicit energy equation, a Newton-Raphson root finder is 
used, supplemented by a binary search algorithm for occasional zones where the tabular 

nature of the cooling curve prevents the Newton-Raphson algorithm from converging. The 
empirical heating function from stellar energy input used in MF99 is turned off in our study, 
Instead, the sharp cut-off in the cooling function at T = 10^ K prevents the background 
atmosphere from spontaneously cooling. 

3.3. Tracer Field and Its Performance 

We use a tracer field c both to follow metal enriched gas and to turn off radiative cooling 
in the hot bubble interior, in order to prevent mass numerically diffused off the dense shell 
from spuriously cooling the interior. The radiative cooling time of the interior is much longer 
than the dynamical time of the bubble for the galaxies of our study, so interior cooling is 
physically unimportant to the dynamics of the bubbles (Mac Low & McCray 1988). We 
set the cooling to be linearly proportional to the value of the tracer field, with the metal- 
enriched, interior gas injected in the source region having c — 0; hereafter we designate this 
as "colored" gas. We do not turn off Compton cooling in the bubble interior because it is 
only linearly dependent on the density. We advect a function 



/(c) = tan[0.997r(c - 0.5)] 



(12) 
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of the tracer field with exactly the same algorithm as the density, rather than the tracer field 
c itself, so that c will show a sharp transition at the interface, even if /(c) becomes quite 
smooth as the result of numerical diffusion (Yabe & Xiao 1993; MF99). 

MF99 used the tracer field successfully to model modern, low-density dwarf galaxies, 
but we find that the bubbles in our high-redshift galaxies are severely poisoned by numerical 
diffusion, despite the sharpening of the edges by the tangent transform. Numerical diffusion 
occurs due to the failure in advection over a large density jump across a small number of 
zones. 

The cooling time of the shells formed by initially adiabatic shocks is tc — 3kTs/4:noA{Ts), 
with postshock electron number density Auq and shocked temperature Tg = {3fis/iQk)Vgf^g^^, 
where Hs is the mass per particle in the ionized ISM, and v shock is shock velocity given by a 
similarity solution of a bubble in a uniform density medium Hq by Weaver et al. (1977). Wc 
use the analytic estimate of the Gaetz and Salpeter (1983) cooling function between 10^ and 
10^ K from Mac Low & McCray (1988), 

A(r) = (1.0 X 10-^^ erg cm^ s'^) (^^) (J-^ . (13) 

This approximation agrees with the primordial cooling function of Sutherland & Dopita 

(1993) that we use for our simulations within a factor of a few in the temperature range. 
Then, we find the instantaneous cooling time of the shocked ambient gas to be very short in 
the early, high-density universe; 

T,, = (0.19 Myr)Ly9 V-^^(Z/0.1Ze)-i oc (1 + z)-'-'' (14) 

for a uniform density medium no oc (1 + z)^, so we expect the formation of dense shells with 
density n oc noAi'^ where Ai is the Mach number. 

We show the results of simple test simulations of the performance of the tracer field in 
Figure 2. Our standard model for this test has fixed background density chosen at z = 13 
Pbg = 200po = 2.5 X 10^^^ g cm~^, which is an average halo gas density at that redshift, 
mechanical luminosity Lm = 10^° erg s~^, and uses the primordial radiative cooling function 
by Sutherland & Dopita (1993). The zone size is 4 pc. We turned off Compton coofing for 
this test because it unnecessarily slows down the computation, but the degree of numerical 
diffusion that we find only becomes worse with the inclusion of Compton cooling. 

With the standard model, the tracer field already fails to pick up 25% of initially marked 
interior gas by t — 5 Myr after the onset of starburst (at time step ~ 5000). The same tracer 
field performs perfectly if the cooling is completely turned off (the adiabatic model: CO), 
while it does worse if the cooling strength is hnearly increased by a factor of 100 (the extra 
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cooling model: ClOO). The CO model does so well because the adiabatic shell is fully resolved 
over ~ 20 cells with density p = 4 x 200 p^g, while the isothermal shells arc unresolved, with 
density p < 200pbgA4^. As the cooling strength increases, the peak density approaches 
its theoretical prediction, and the numerical diffusion gets worse. In the low-resolution 
simulation LR with 12 pc zones, the performance of the tracer field is even worse due to 
the failure to resolve the shell. Numerical diffusion increases in our models as unresolved, 
high-density, swept-up shells fragment after the bubbles blow out of the disks, increasing the 
area of contact between high and low density gas. 

Wc find that the degree of numerical diffusion near unresolved shells docs not improve 
substantially by changing the advection scheme from second-order, van Leer method to 
third-order Piecewise Parabolic Advection (PPA) method (Clarke 1988; Colella & Woodward 
1984). We also find that changes in the size of source regions do not influence the performance 
of the tracer held. 

The tracer fleld controls cooling, preventing it in gas marked as interior ejecta. When 
hot gas is mixed into regions colored as external, it is allowed to cool. As mass from the shell 
has usually numerically diffused into this gas, it cools unphysically quickly, so as diffusion 
proceeds, the interior energy of the bubble is artiflcially dissipated. Our goal in this study is 

to compute metal/mass ejection and energy transport efficiencies by simulating the bubble 
dynamics in high-rcdshift galaxies with the tracer field. However, our simulations only 
yield strict lower limits on the efficiency of metal, mass, and energy feedback because of 
the numerical diffusion. To solve this problem we will eventually need cither high enough 
resolution to cleanly resolve shells, perhaps using adaptive mesh refinement, or else use of a 
true contact tracer as was done by Mac Low et al. (1989). 



4. Galaxy and Star Formation Parameters 
4.1. Galaxies 

We model objects representing 1-3(7 peaks in the ACDM density perturbation spectrum 
(see Fig. 3). We focus on second generation systems with virial temperature Ty > 10^ K 
that can cool through line cooling by H and He. We show the galaxies we model with halo 
masses Mh = 5 x 10^-5 x 10^° M© that are expected to collapse at z = 3, 5, 8, and 13 in 
Figure 3. We make a list of the properties of our model galaxies in Table 1. 

In these second generation objects, the cooling time is much shorter than the dynamical 
time, so the gas will cool nearly isothermally at ~ 10^ K and rapidly ffow to the centers. 
We assume that it settles in tightly, bound, rotationally supported disks (Navarro & White 
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1994) without much fragmentation (Kashhnsky & Rees 1983; Oh & Haiman 2002). Such 
disks are mostly gravitationally stable, but a large fraction of out-of-equilibrium, residual, 
electrons expected during initial line cooling drives H2 formation through gas-phase reactions 
catalyzed by H~ (Mac Low & Shull 1986; Shapiro & Kang 1987). Radiation from the H2 
allows the gas to cool further and collapse (Oh & Haiman 2002). 

Figure 3 also shows the Jeans mass for gas in photoionization equilibrium with the 
background UV radiation with T ^ 2 x 10"^ K. We chose two galaxies with = 5 x 10^ Mq 
and fc ~ 30 — 40 km s~^ at z = 8 and 13, assuming the effect of the UV background 
radiation to be negligible at ^ > 6. The cooling of gas in halos with Vc < 30 km s~^ is 
suppressed in the reheated IGM (e.g. Bond, Szalay, & Silk 1988; Mciksin 1994; Efstathiou 
1992; Thoul & Weinberg 1996). Therefore, we would not expect starbursts in such small 
halos at 2; = 8, if they are found in the IGM fully ionized by that time. The results in these 
halos should be interpreted in that hght. We also chose three galaxies with Mh = 5 x 10^ M© 
and ~ 40 — 60 km s~^ at ^ = 3, 5, and 8. These galaxies fall below the Jeans mass, but 
hydrodynamic studies of dwarf halo formation show that ~ 50% of baryons within the virial 
radius can still cool as long as > 50 km s"^ (Vedel et al. 1994; Thoul & Weinberg 1996; 
Navarro & Steinmetz 1997). We can still expect moderate starbursts in these halos. 

As a comparison, we also plot in Figure 3 the minimum halo mass for first generation, 
or Population III, systems with T^, < 10"^ K that can only cool through molecular hydrogen 
formation, as computed by Tcgmark et al. (1997) for a standard CDM model. We do 
not select any of these first generation objects because their star formation efficiency is 
expected to be very low due to the presence of internal and external photo-dissociating 
radiation (Omukai & Nishi 1999; Glover & Brand 2001; Haiman et al. 2000). We are far 
from realistically computing the efficiency of cooling and star formation even in second 
generation systems that we study, and cannot definitively predict in which halos starbursts 
will occur at a given redshift. Instead, we perform a parameter study that explores a large 
range of dwarf galaxies to understand the physical processes of stellar feedback in different 
environments. 

In Figure 1, we plot our chosen galaxies on top of the pressure confinement criterion 
given by equation (8) for different values of i/f^. With = 0.1, equation (8) predicts that 
any bubble in a dwarf galaxy will escape to the IGM unimpeded by the ram pressure of 
infalling gas, unless the bubble has lost so much of its energy that z/ < 0.1 and the galaxy 
formed from an extremely large density perturbation (see Fig. 3). With = 0.001, most 
second generation systems with T > 5 x 10^ K are predicted to be confined by ram pressure 
with u < 0.5. The most interesting case in Figure 1 is with /* = 0.01. If the fraction of energy 
preserved in the bubbles u < 0.2, the bubbles in halos with Mh > 5 x 10^*^ Mq from ~ 2a 
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peaks at z — 3-5 will be confined by the ram pressure of the infall. If < 0.1, even bubbles in 
smaller halos with ~ 5 x 10^ Mq will begin to be confined by the ram pressure. However, 
we can expect the bubbles in halos with A4 < 10^ (T < 5 x 10'' K) to escape to the 
IGM, unless most of their energy is lost to the surroundings. Figure 1 suggests high metal, 
mass, and energy feedback efficiencies in small high-redshift dwarf galaxies with /* = 0.01 
and larger galaxies with /* 3> 0.01. In §5, we compare the results of our simulations with 
these predictions, and try to constrain i/. 



4.2. Disks 

In our simulations, we set up a disk with a uniform surface density profile, E = M^/vri?^, 
where R^, is the size of the disk, as in MF99. We choose a uniform surface density profile 
for numerical simplicity, although an exponential surface density distribution is observed in 
local blue compact dwarf galaxies (van Zee et al. 1998). However, bubble evolution is mainly 
governed by the vertical distribution of gas within the central region of the disk, where a 
constant surface density is a good approximation. 

To compute the size of a galaxy as a function of halo mass and collapse redshift we 
follow Mo, Mao, & White (1998), using an isothermal halo profile. The size of a disk with a 
uniform surface density is 

H. ^ I.(^)XR.. (15) 
where is the virial radius of a halo, 

R^ = (16 kpc)(l + z)-' [y^) i^ohT'^', (16) 

where we use the same approximation for the ACDM model as in §3.1, and we set the 
dimensionless spin parameter A — 0.05, approximately at the peak of its distribution (e.g. 
Cole & Lacey 1996). We assume that a fraction of halo mass, rud — Md/M^, settles into a 
non-self gravitating, rotationally supported disk in a singular isothermal halo with a fraction 
of halo angular momentum, — Jd/ Jh- This method provides a good fit to the observed 
size distribution of galactic disks for jd/ma — 1 (Mo, Mao, & White 1998). 

The baryonic content of each halo is Mg — {^lb/Qo)Mh^ the amount of gas within the 
virial radius R^ if the gas simply follows the dark matter. We choose the mass which goes 
into the disk Md by subtracting the result of the infall solution described in §3.1 from 
Mg. Typically, ~ 0.5Mg, but it varies by 10% among models with different cooling 
histories. The disks containing ~ 50% of cooled baryons are appropriate for dwarf galaxies 
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with Vc ^ 50 km s~^, but not for those with Vc ~ 30 km s~^, if the infall gas was initially 
reheated to a few xlO^ K (see §4.1). In halos with Vc ~ 30 km s~^ a,t z < Zion, most of 
the halo is pressure-supported, and only small fractions of the gas mass Mg can cool to 
form disk-like objects. The evolution of bubbles in such galaxies will only differ within the 
virial radius from our standard models with disk masses as high as Md ~ 0.5Mg. The 
qualitative results on metal, mass, and energy feedback efficiencies should not differ much, 
as long as disks or disk-hke objects with stratified ISM are assumed to form, because the 
bubbles blow out of them quickly to interact with the halo gas and the IGM. However, star 
formation efficiency /* should be strictly < 0.01 in the halos with Vc ~ 30 km s~^ if the 
collapse redshift z < Zion- 

We initially set up the disk in hydrostatic equilibrium with a potential field due to the 
sum of the dark matter halo at the collapse redshift and the potential of a cold, thin disk 
representing earlier stars formed in fragmented gas or in the disk prior to any star burst. The 
dark matter potential computed by the infall code evolves as a function of time, but we do 
not modify the initial potential within in our two-dimensional models. Since our dark 
matter profile is not realistic in the center, as discussed in §3.1, we must include the cold, 
thin, disk potential to model a disk with an exponential atmosphere that resembles the gas 
distribution observed in the Galaxy and nearby galaxies (Lockman et al. 1986; Neeser et al. 
2002). To do so, we assume the presence of rather massive, thin disks containing 0.03Mg. 
This approximation of a fixed, cold disk potential is probably reasonable, as the dark matter 
still dominates the potential everywhere. 

We assume that the gas in the disk is in the warm, neutral phase with some turbulent 
contribution provided by winds and supernovae from the same small-scale star formation 
that formed the cold, thin disk, yielding an effective sound speed of Cg ~ 10 km s~^ with 
corresponding temperature, T ^ 10^-^ K. Modeling this diffuse, warm medium is important 
for the dynamics of bubbles, since the bubbles will immediately blow out of the cold gas 
disks where star formation is actually expected to occur (Mac Low & McCray 1988). 

4.3. Starbursts 

At the beginning of our simulations, we set up an instantaneous starburst at the center 
of each disk. The expansion of a starburst bubble can be considered to be an extension of 
superbubble dynamics in an OB association in a disk galaxy. The stellar winds create a hot, 
low-density cavity in the ISM, and repeated supernovae excavate a larger hole as they sweep 
the gas into a thin, dense shell (Tomisaka & Ikeuchi 1986, Mac Low & McCray 1988, Mac 
Low et al. 1989). Stellar winds are weaker in low-metallicity stars (e.g. Kudritzki 2002), but 
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the initial few SNe will still create a similar hot cavity in high-redshift galaxies. Once the 
hot cavity forms, discrete supernova explosions generate blast waves that become subsonic 
in the hot interior, and hence can be treated as a continuous mechanical luminosity in 
the study of bubble dynamics (Mac Low & McCray 1988). 

It is, of course, an approximation to assume that all the mass available for the starburst 
turns into stars instantaneously and that the starburst is concentrated in a single spot at the 
center of a disk. Local observations of dwarf starburst galaxies show that multiple starburst 
clumps are scattered around the disks instead (e.g. Vacca 1996; Martin 1998). High-redshift 
disks are much smaller than present-day dwarf galaxies, roughly oc {1 + z)^^, but the 
smallest galaxy that we model extends over a few hundred pc, which is an order of magnitude 
larger than a typical molecular cloud (~ 10 pc). Therefore, we expect the same or at least 
similar starbursts in high-redshift galaxies, too, as long as coohng is efficient. The starburst 
clumps are made up of OB associations or clusters. Since the largest OB associations or 
clusters have masses of a few million solar masses, we could actually think of our starbursts 
as being powered by a single star cluster as long as we do not assume too high a star formation 
efficiency in our model. 

In our study, we use uniform mechanical luminosities computed using the Starburst 99 
model with Z = 0.001 (Leitherer et al. 1999). The Starburst 99 model is based on a power 
law initial mass function which approximates the classical Salpctcr (1955) function, with 
exponent a = 2.35 between low-mass and high-mass cutoff masses of 1 Mq and 100 Mq 
respectively. The average mechanical luminosity predicted by the Starburst 99 model is 
Lm = 2.7 X 10*^° erg s^^ when 10^ M© of gas is instantaneously converted to stars of mass 
between 1 Mq and 100 Mq. Note that approximately the same mass of stars below 1 Mq 
will form at the same time. The Starburst 99 model is based on observations of star-forming 
and starburst regions in the local universe. However, various low to intermediate-redshift 
observations suggest that the initial mass function may have been top heavy in the early 
universe (see the summary by Larson 1998), a suggestion supported by recent studies of 
Population HI stars (Abel et al. 2000; 2002; Bromm, Coppi, & Larson 1999; 2002; Nakamura 
& Umemura 1999; 2001). If the initial mass function is indeed top-heavy at high redshift, 
the fraction of O and B stars which contribute to the mechanical luminosity increases, so the 
conversion of the same amount of gas to stars can provide a higher mechanical luminosity 
than predicted by the Starburst 99 model. Within the uncertainties, we assume that the 
mechanical luminosity of a burst is proportional to an amount of gas f*Mg available for a 
starburst as follows: L„ = 2.7 x 10^V*(^g/10^ Mq) erg s-^ 

Our mechanical luminosities range from 10^^ < < 10^^ erg s~^ produced when stars 
form from 6-7x10^ Mq of gas in the smallest halo with = 0.001 and 6-7x10^ Mq in the 



-20- 



largest halo with = 0.1. These mechanical luminosities are equivalent to a supernova with 
energy of 10^^ erg exploding every ~ 17, 000 yr to every ~ 1.7 yr, respectively. Energy input 
is continued for a period of 50 Myr, roughly the lifetime of the smallest B star to go supernova 
(McCray & Kafatos 1987). Most cosmological simulations with supernova feedback employ 
> 0.1 to test its maximum effects. 

We consider starbiu^st efficiencies in our galaxies of 0.001, 0.01, and 0.1. Numerical 
simulations of first star formation suggest ^ 0.001 in metal-free primordial clouds (Abel 
et al. 2000; 2002; Bromm, Coppi, & Larson 1999; 2002). Oh & Haiman (2002) estimate an 
initial star formation efficiency by molecular hydrogen cooling alone to still be = 0.001, 
even in second generation disks. The star formation efficiency in such disks is not well 
constrained because of the presence of internal, photo-dissociating, UV radiation, but we 
can expect metal cooling to enhance the star formation efficiency, if the halo gas is pre- 
enriched at least to Z — 5 x 10~^ Zq by previous generations of stars (Bromm et al. 2001; 
Oh & Haiman 2002; Wada & Venkatesan 2003). Therefore, the star formation could actually 
be rather efficient in second generation objects. We think that /* ~ 0.01 may be a reasonable 
guess for the initial major starburst in a galaxy. Starburst efficiencies as high as /* = 0.1 are 
probably not realistic. Since we are performing a parameter study, however, we do consider 
them, in order to clearly delineate the starburst efficiencies with which efficient metal/mass 
ejections and energy transport are expected in a halo at a given redshift. 



5. Results 

We have explored the effects of repeated supernovae with mechanical luminosity linearly 
proportional to the amount of mass converted to stars, /*M^, in seven examples of early 
dwarf galaxies. Figure 4 shows time histories of the density distributions of our model with 
Mh — 5x 10^ Mq ai z = 8, with starburst efficiencies /* ranging from 0.001-0.1. The bubbles 
shown in Figure 4 are representative of our other models in a number of ways as we now 
describe. (Figures 5a-c show the density distributions of all of our models a,t t — 100 Myr.) 

The bubbles quickly begin to accelerate vertically through the stratified ISM of the 
galaxies, and blow out of the disks by i ~ 10 Myr after the onset of the starbursts. With 
sufficiently high central mechanical luminosity (/* = 0.01-0.1), the bubbles sweep up the 
entire disks by i = 30 Myr. This behavior occurs in all of our galaxies, because our high- 
redshift disks are very small, less than ~ 1 kpc in size, and are highly stratified in their deep 
potential wells. Most of the ISM remains bound, so these disks are not yet "blown-away" in 
the sense used by De Young & Heckman (1994) and MF99. In fact, we will show below that 
most of the ISM remains bound to the halo potentials if < 0.01. 
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The internal termination shocks outside the freely expanding winds and ejecta are clearly 
seen at t = 10-30 Myr, surrounded by hot, pressurized regions of shocked winds and super- 
nova ejecta. The hot, pressurized regions drive shocks into the surrounding ISM and halo 
gas, sweeping it up. In these high-redshift, dense galaxies, gas behind the outer shocks tends 
to cool efficiently, so that the shocks are usually isothermal. Outside the galactic disks, the 
shocks are no longer isothermal. Figure 6 shows the formation of dense shells, cooled from 
Tshock ~ 10^ K, driven by the hot, pressurized gas at T f« 10^ K behind, as the bubble blows 
out of the disk at t = 10 Myr, with = 0.01 (middle panel in Figure 4). The swept- up gas 
forms dense, cold, thin shells except in the cases with /* = 0.1, when the shocks are so fast 
that cooling remains inefficient. 

Shell formation is extremely sensitive to the resolution employed, and is not fully re- 
solved in our simulations (see Figure 6 and also the resolution study of shell formation in 
Fujita et al. 2003). Resolving the shell is not important to following the overall dynamical 
evolution of bubbles driven by the thermal energy of the hot pressurized regions (Castor, 
McCray, & Weaver 1975; Weaver et al. 1977), but is important to follow the details of 
shell fragmentation due to Rayleigh- Taylor instability when the bubbles accelerate. Figure 4 
shows examples of shell fragmentation. The development of the Rayleigh- Taylor instability 
in our model is limited not only by the resolution, but also by the assumption of azimuthal 
symmetry. Mac Low et al. (1989) pointed out that in an axisymmetric blowout, the typical 
spike and bubble structure of the Rayleigh- Taylor instabihty is hmited to rings. Thus the 
detailed structure of the fragments will be different in 3D. However, the instabilities seen in 
our simulations effectively allow hot gas to accelerate beyond the cold, swept-up shells, which 
is the qualitative result of interest. Figure 4 shows that for /* > 0.01 the hot, metal-enriched 
gas expands into the IGM well beyond the virial radius of the halo, leaving the fragmented 
shells behind, which eventually fall back into the center. The hot gas is confined in the 
vicinity of the disk if = 0.001. Similar behavior can be seen in Figure 5a-c for all of our 
model galaxies. 

We now discuss our quantitative results. We compute metal and gas ejection and energy 
transport efficiencies in §5.1; discuss the effects of external ram pressure from infalling gas 
on the evolution of the bubbles and so on the feedback efficiencies in §5.2; and estimate in 
§5.3 how long it takes for swept-up but still bound mass to fall back to the center of a halo, 
where it will be available to fuel the next star burst. 
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5.1. Feedback Parameters 

In order to quantify feedback, we measure metal ejection efficiency .^met, mass ejection 
efficiency ^, and energy transport efficiency (. In each case, we measure the content of gas 
that has been accelerated to outward velocities higher than the local escape velocity defined 
by the gravitational potential. Because of the limitations of our model, we actually measure 
the efficiencies when the bubbles reach 2Ry, where we assume the low-density IGM to begin. 
This choice is rather arbitrary, and is only motivated by the idea that the bubble should 
be well outside the halo virial radius Ry, in order to see the effect of the ram pressure of 
inf ailing gas. 

Recall that we use a one-dimensional, spherical background potential (§3.1). Simulations 
show filamentary structures are anchored in spheroidal halos (Zhang et al. 1998). Thus the 
spheroidal approximation is a good one within a virial radius, but becomes progressively 
worse beyond. In reality, bubbles such as those we model would possibly find their way out 
to low-density regions between the filaments. When the bubble finds such a funnel for escape, 
we expect that any gas carrying supernova energy and metals that has velocity greater than 
the escape velocity of the halo potential will escape nearly unimpeded into the IGM. 



5.1.1. Metal Ejection Efficiency 
The metal ejection efficiency was defined by MF99 as 

^met,M^' ~ ^c,esc/MsN, (17) 

where M^^esc is the mass of metal-enriched gas moving with speed greater than the local 
escape velocity and colored by tracer field with c < 0.6; and Msn is the total mass injected 
at the center with c = initially, which accounts not only for winds and supernova ejecta, 
but also for the mass evaporated off the swept-up shell by thermal conduction. The colored 
mass Mc traced is not sensitive to the choice of the cut-off number (0.6) in the absence of 
numerical diffusion. The mass Msn is found inside the hot, pressurized, region in which 
metals are assumed to be mixed well. Equation (17) is based on the assumption that the 
injected mass Msn is preserved at the time of its measurement, without being mixed into 
the ambient medium physically or numerically. However, the numerical diffusion severely 
poisons the colored gas as discussed in §3.3, and our tracer field only picks up ~ 60% of 
originally colored gas within 5-10 Myr after the onset of the starbursts. 

At the same time, we observe a significant amount of physical mixing between the 
metal-enriched, colored gas and the ambient gas. We show in Figure 7 an example of metal 
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mixing inside a bubble in a halo with Mh — 5x 10^ Mq at 2; = 8. The dense, swept-up shell 

of ambient gas fragments due to Rayleigh- Taylor instability when the bubble blows out at 
t ^ 10 Myr. Figure 7 shows the density distribution of the bubble at t = 19 Myr, and shows 
that the fragmented shells do mix with the hot interior gas very effectively. This physical 
mixing of gas accelerates numerical diffusion further. Unfortunately, we have no systematic 
way of tracking physical metal mixing, nor of distinguishing it from numerical diffusion, 
but it is important to quantify the fraction of the escaping gas enriched with metals to any 
extent. 

Therefore, it is not appropriate for us to use equation (17), since we severely under- 
estimate the metal ejection efficiencies when we can not track down the total amount of 
originally colored mass. Instead we define the metal ejection efficiency as 

Uet = M,^,JM,, (18) 

where now the color cut-off, c < 1 — e, with e chosen small. In this way, we are sure to 
pick up any metal-enriched gas (that is, gas that was injected with c = 0), and still obtain 
a strong lower limit to the actual metal ejection rate, since Mc ^ Msn- Figure 8 shows 
the colored mass as a function of e for the halo with M^, = 5 x 10® at 2; = 8 for 
star formation efficiencies = 0.001, 0.01, and 0.1. The masses Mc were measured when 
the bubbles reached 2R^ at t = 100 Myr when /* = 0.001 and 0.01 and at i = 50 Myr 
when /* = 0.1. The mixing between the colored gas and the ambient gas proceeded far by 
t — 100 Myr for the bubbles with = 0.001 and 0.01 (see Figure 5a), so the values of Mc 
are insensitive to the choices of e as long as they are very small. The asympotic value of Mc 
as e decreases is not reached for the case with = 0.1 because it was measured at an earlier 
time. However, the bubble is so powerful as to reach 2i?^ by that time, so any values of Mc 
will yield irnet ~ 1 with 10"*^ < e < 0.6. We choose e pa lO'MO-^ to enforce < M^, the 
total mass of baryons in the halo. 

The metal ejection efficiencies ^met are nearly zero in all the halos that we study when the 
star formation efficiency /* = 0.001 and close to 1 with = 0.1 (Table 2). These results are 
consistent with our prediction in §2 that in the presence of infall ram pressure, the bubbles 
in second generation objects can escape to the IGM if = 0.1, but will be confined within 
the halo potentials if /* = 0.001 and v < 0.5. We find that the metal ejection efficiencies 
are largely determined by the value of /* that we assume. Recall that /* = 0.1 is likely to 
be unrealistic for instantaneous starburst in second generation systems. 

The most interesting results are therefore the values of ^rnet with = 0.01. Metal 
ejection efficiencies are low in halos with Mh — 5 x 10^° Mq at 2; = 3 and 5 and even with 
Mh — 5x 10^ M0 at 2; = 8. However, high metal ejection efficiencies ^met ^ 0.5 occur in halos 
with virial temperatures T^, < 1 x 10^ K or circular velocities Vc < 50 km s~^ at all redshifts 
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we studied. These galaxies may make significant contributions to the metal-enrichment of 
the IGM. On the other hand, larger galaxies with Vc ~ 100 km s~^ allow efficient cooling, so 
star formation efficiencies 3> 0.01 may be plausible. These galaxies could then produce 
and expel more metals to the IGM. 

If we assume that each supernova produces ~ 3 of heavy elements (Ferrara & Tolstoy 
2000), a halo with dynamical mass can expel metals with mass 

Mrr^et ~ (8 X 10^ Mq) ( ( ( ^ ^ . (19) 

If we further assume that those metals mix uniformly, we can estimate the average metallicity 
of the universe from halos at a given redshift z; 

^ " ^ iM..,) I lMpc-3 J J [r^) IsTTo^ j ' 

(20) 

where Nh{Mfi, z)dMh is the comoving number density of halos at z with masses between 
M/j and + dM^, predicted from Press- Schechter theory (Press & Schechter 1974, Bond 
et al. 1991) in our ACDM model, and M^ax and Mmin{z) are the maximum and mini- 
mum halo masses to be incorporated into the calculation. We set M^nmiz) — 5.8 x 10^ 
Mq[{1 + z)/9]~^^'^, corresponding to Vc — 30 km s~^. Although the metals may well be 
distributed highly inhomogeneously, equation (20) provides a lower limit to the average 
metallicity of intergalactic clouds where the expelled metals can reach. 

Assuming ^metf* to be constant, we plot the average metallicity Z oc Cmetf* (eq. 20) as 
a function of mass Mmax at 2; = 3, 5, 8, and 13 in Figure 9. We chose ^metf* = 0.005, which 
means for example, = 0.01 and C,met = 0.5 for halos with Vc ~ 50 km s~^ and = 0.05 
and ^met = 0.1 for halos with Vc ~ 100 km s~^. Although we neglect the dependence of ^met 
on halo properties for a given Figure 9 shows a qualitative trend of metal contributions 
from halos at different redshifts. 

In Figure 9, the metal contribution seems to be poor from halos at 2; = 13, but rather 
substantial from halos at z = 8. The average metallicity at z = 8 from halos with v,. < 50 
km s~^ alone is already approaching Z ^ 10"^ Zq, the lower limit of the metallicity observed 
in low-density Lyman a clouds at 2; ~ 3-5.5. Such halos correspond to ~ 2a peaks in our 
ACDM primordial density perturbations. The early enrichment at 2; ~ 8 may explain the 
absence of kinematic disturbance observed in low-density Lyman a clouds at 2; ~ 3 (Ranch 
et al. 2001b; Madau, Ferrara, & Rees 2001), by which time the kinematic effects of the 
wind may have dissipated and the gas resettled onto the filaments. Note however, the metal 
contributions from halos Vc ^ 50 km s~^ are a factor of a few larger at lower redshifts, since 
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they are more abundant (1-2(7 peaks) with larger halo masses, and so can produce more 
metals with the same ^metf*- Figure 9 also shows that the larger galaxies with Vc ~ 100 km 
s~^ begin to make significant contributions, yielding Z > 10~^ at ^ ~ 3-5. The high Hi 
column density Lya clouds (logA^j^j > 14.4) with Z ^ 10~^ Zq (Cowie et al. 1995) might 
have been recently enriched by such dwarf galaxies then. This picture of recent enrichment 
is also supported by large turbulent motions {vrms ~ 70 km s "^) observed in high column 
density CiV systems (Ranch et al. 2001b). 

Our estimates in Figure 9 based on our numerical simulations suggest that both dwarf 
starburst galaxies with Vc ~ 30-50 km s~^ and with v,. ^ 100 km can be very cfiicient 
in spreading metals, but to difi^crcnt regions of the IGM (low and high density clouds) at 
different times (z ~ 8 and ^; ^ 5). We find that only galaxies collapsing from 1-2(7 peaks 
are abundant enough to make significant contributions to the metal enrichment of the IGM. 

Metal ejection efficiencies ^m,MF ~ 1 were found by MF99 in all of the local dwarf 
galaxies that they modeled, even with star formation efficiency = 0.001, while we found 
neghgible metal feedback with = 0.001 in high-redshift galaxies. We showed in §5.2 that 
the difference between our results and those of MF99 is not the higher densities p oc {l-\-zf, 
because bubbles quickly blow out of the smaller disks with smaller scale heights expected 
at higher redshift. Rather, it is the ram pressure from cosmological infall that makes the 
difference. 



5.1.2. Gas Mass Ejection Efficiency 
We define the efficiency of gas mass ejection as 

^ (fi,M)M,' ^ ^ 

where M^sc is the mass moving with velocity v > Vesc, and the denominator is the mass of 
baryons initially within Ry, as we assumed they initially simply follow the dark matter in 
our models. We fist ^ for each of our models in Table 2. 

The mass Mc,esc of escaping, metal-enriched gas is a negligible fraction of Mesc- The 
bulk of the mass escaping is some combination of ISM and infalling halo gas. Since we 
do not systematically track the disk gas in our simulations, we identify gas with angular 
momentum as disk gas, although there is some transfer of angular momentum between the 
disk and ambient gas due to numerical diffusion and physical mixing. Our accuracy is 
fimited, so we can only observe whether > 50% of the ISM is moving faster than the escape 
velocity or not. Halos in which this occurs are marked as having had their ISM blown away 
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in Table 2. 

Table 2 shows that the mass ejection efficiencies ^ >l with = 0.1 and are nearly zero 
with = 0.001, consistent with the metal ejection efficiencies. More than half of the disk 
gas escapes the potentials when /* = 0.1 in all the halos except one. Blow away is possible 
with /* = 0.1 because high-redshift galaxies are very small and easily swept away by powerful 
bubbles: however, om result should be interpreted carefully because of our assumption of 
an instantaneous starburst at a single site (see the discussion in §4.3). With = 0.01, 
the mass ejection efficiencies ^ can be large for the same parameters when C,met is large, but 
most of the escaping mass is from the halo, while the ISM remains bound. (This agrees 
with MF99 and D"Ercole & Brighenti [1999], who found that ISM loss is inefficient so long 
as it is distributed non-spherically.) Thus, a substantial fraction of infalling halo gas can 
acquire enough energy from the bubbles to escape the potentials. For example, halo gas 
mass comparable to {Qb/Qo)Mh escapes from the halo potentials with Mh — 5 x 10^ Mq 
at 2; = 8 and Mh = 5 x 10^ M0 at z = 3, even with /=„ = 0.01. These results have two 
interesting implications for galaxy formation: further cooling in the halos can be suppressed 
to a certain degree, which we discuss in §5.3; and the turned-around halo gas will carry the 
energy of supernova explosions to other sites of halo formation, perhaps preventing their 
collapse, which we discuss next. 



5.1.3. Energy Transport Efficiencies 



In order to compute the energy transport efficiency we add up the kinetic energy K 
and thermal energy U carried by any gas with v > v^sc^ subtracting the potential energy 
required to escape the local potential well: 



c 



i^e.c + C/e.c-(^>(2i?.)-$(r)) 



Lmtf 



(22) 



where $(r) is the potential of a gas parcel located at r, $(2i?^) is the potential at 2i?^ 
ft-om the center, and the input energy is L^tf. We exclude thermal energy carried by any 
swept- up shell of ISM material and infalling gas, if the Hubble time at a given redshift tniz) 
exceeds the coohng time of a shell with radius Rg, 



Tsh = (78000 yr) 



l + z 



0.01 



1.7 



Rv{Mh,z) 



0.3 



z 



(23) 



where we assume that the density profile is that of an isothermal halo [Eq. (6)], the bubble 
is driven by the residual energy i^L^tf with star formation efficiency and Ry is the virial 
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radius for a given halo Mh at a given redshift z. We used equation (13) for the coohng rate 
A(r). In fact Tsh < tniz) unless = 0.1 at z < 8. 

We can not systematically distinguish between energy carried by hot, metal-enriched gas 
and swept-up, ambient gas because of numerical diffusion. However, we can safely compute 
the amount of thermal energy carried by escaping shells of swept-up gas with a conservative 
color cut-off, c > 0.6, because the thermal energy is carried by the gas which is recently 
shock heated and therefore not mixed with ejecta by numerical diffusion. In Figure 10, we 
plot the fractions of energy carried by ejecta and by ambient gas with c > 0.6 in a halo 
with Mh = 5 X 10^ Mq at z = 3. It shows that swept-up halo gas carries most of the 
energy, in kinetic form. Even if we raise the color cut-off to c > 1 — e, the kinetic energy 
of the ambient gas remains much larger than that of the interior gas. After the bubble 
blows out at i 10 Myr, the thermal energy of the hot interior is effectively transferred to 
the surrounding gas. The bubble reaches 2Ry by i ~ 100 Myr, and the fraction of energy 
escaping C 0.18. Energy transport efficiencies ( for all the models are listed in Table 2. 

We note that numerical diffusion could accelerate the coohng of the bubbles. However, 
a resolution study of energy in models with = 0.01 and — 5 x 10^ Mq a,t z — S using 
standard and doubled resolution shows no more than 30% differences in energy (Fig. 11). The 
largest differences are in shell thermal energy, which is least affected by numerical diffusion. 
This suggests that the dissipation of bubble energy seen in our simulations is physical, not 
dominated by numerical diffusion. 

The transfer of the blast- wave energy into the IGM is effective, with C — > 1, in models 
with /* = 0.1 and negligible if = 0.001, as expected from our discussion of the metal and 
mass ejection efficiencies. The energy transport efficiencies ( are roughly proportional to the 
mass ejection efficiencies ^, because the escaping mass is carrying most of the energy. Both 
seem to be largely determined by With = 0.01, the energy transport efficiencies are 
C ~ 0.1-0.3 in halos with Vc ^ 50 km s~^. Larger halos with Vc ~ 100 km s^^ have more 
efficient cooling, and likely therefore have higher star formation efficiencies /* ^ 0.01, so 
that they can transport more total energy into the IGM. 

This kinetic energy feedback could stir intergalactic gas and support it against gravi- 
tational collapse in other halos. Recent numerical work has demonstrated that supersonic 
turbulence can delay gravitational collapse for many free-fall times in turbulent molecular 
clouds (Gammie & Ostriker 1996; Klessen, Heitsch & Mac Low 2000). However, turbulent 
energy decays quickly, in a time roughly proportional to the crossing time tcr — L/vrms of 
a turbulent region of size L with flow having rms velocity Vrms (Mac Low 1999; Elmegreen 
2000). The distance reached by the turbulent motions is only ~ L-\/2 after half the energy 
is lost at tcr (Avila- Reese & Vazquez-Semadeni 2001). 
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We now crudely estimate whether turbulence driven by starburst bubbles in the very 
first dwarf galaxies could also delay and suppress the formation of subsequent dwarf galaxies 
with the same or larger mass Mg > Mj. We can estimate how much turbulent energy 
a single dwarf starburst galaxy can supply to surrounding overdense regions with baryonic 
density p, by dividing the escaping energy by the volume of the bubble at Ri,: 



2 < (I,„f/ > 
(4/3)7rpi?g 



(40 km s"^) 



1/2 



0.1 



0.01 
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(24) 
(25) 



where we used equation (7) with tf — 50 Myr. This depends on energy ejection efficiency 
star formation efficiency /*, and halo virial radius R^ — Ry{Mh, z) from equation (16). Note 
that the average density of a gravitationally bound object at turn-around p/p — 97r^/16 ~ 5.6 
from the spherical collapse model. Based on our simulations, we estimate v^ms ~ 40 km s~^ 
with = 0.01 from galaxies with < 50 km s~^, assuming p/p = b and = 5i?.u. The 
same level of turbulence could result from efficient star formation ~ 0.05 yielding C, ~ 0.02 
in galaxies with ~ 100 km s^^. Recall that the kinetic energy of the wind is carried by 
the swcpt-up halo and intcrgalactic gas and is transformed to turbulence through dynamical 
instability. Therefore, equation (25) underestimates the amount of turbulence at the edge 
of the bubble, but overestimates it within the bubble. 

In order to delay and suppress cooling in halos with fc, the turbulence with Vrms ~ Vc 
must be driven and maintained by repeated starbursts occurring on a dissipation timescale 
of the turbulence in equation (25): 



^cr Rvir / ^rms ^ 0.3Gyr 



< tH{z) ^ 1.9 Gyr 



4 



-3/2 



-3/2 



Note that some fraction of cooling in the halos with Vc can be suppressed by Vr, 
since the turbulent pressure delays the turn-around of gas from the Hubble flow. 
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Vc/2, 



The Press-Schechter mass function for our ACDM model predicts the mean separation 
of dwarf galaxies corresponding to 2a peaks to be ~ 10i?„, independent of mass and redshift, 
when they are assumed to reside in overdense filaments with p/p ~ 5-10 (corresponding to 
a volume fiUing factor of a few percent). This is hkely to be an overestimate, since 2a peaks 
will be clustered. In addition, the mean separation of all the dwarf galaxies at a given epoch 
will be smaller, so we assume bubbles from these galaxies collide with other bubbles at ~ 5i?^ 



-29- 



on average. Starburst bubbles are powerful enough to reach distances much greater than a 
few virial radii hj t = 50-100 Myr in our models with < 50 km and = 0.01 and in 
all models with = 0.1. Therefore, we can only expect these sorts of halos, which typically 
form from l-2o" peaks, to provide turbulence of the order ~ 40 km s^^ to the IGM. Such 
starburst-driven turbulence may delay and suppress the cooling of the gas in future halo 
formation sites. 



The turbulent Jeans mass (Chandrasekhar 1956) is 



Mjt ^ 2.4 X lO^^Me 
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, 16 km s ^ / 3 \ 16 km s ^ 

where = 16 km s~^ is the sound speed of the IGM with T = 2 x 10^ K. The Jeans mass 
is 5.4 times larger than the thermal Jeans mass alone for halos with Vrms — 40 km s~^, 
corresponding to 120 km s". 

This mechanism is analogous to the suppression of dwarf galaxies by UV background 
radiation (Bond, Silk & Szalay 1988; Meiksin 1994; Vedel et al. 1994; Thoul & Weinberg 1996; 
Navarro & Steinmetz 1997; Iliev, Raga & Shapiro 2003), except that the gas is energized 
by turbulent flows rather than thermal pressure due to ionizing radiation. Semi-analytical 
studies of galaxy formation suggest that suppressing the fraction of gas that can cool and 
form stars below ~ 10-20% in halos with Vc ~ 50-100 km s~^ can reproduce the shape of 
the observed galaxy luminosity function (e.g. Cole ct al. 1994). 

Some cosmological simulations have demonstrated that the galactic outflows that they 
model, though not based on superbubblc dynamics, can suppress the formation of dwarf 
galaxies by heating up the gas by shocks (Thcuns ct al. 2002) and by stripping the gas from 
the potentials by their momentum ("baryonic stripping" Scannapieco et al. 2001). We agree 
that starburst-driven turbulent flows may play an important role in regulating the formation 
of galaxies, and so may be a key to the long-standing problems of overcooling and angular 
momentum losses in a CDM universe. It appears useful to model such starburst-driven 
turbulence in more detail in cosmological simulations, in order to address such questions. 

In terms of the locally induced metallicity Z, the rms turbulent velocity is 

1/2 / ^ ^ 1/2 
Imety \10~^ Zq 



(60kms-)(^) [j^^] . (28) 



The trend toward larger galaxy masses dominating metal ejection at later times (see Fig- 
ure 9) suggests that Vrms may increase with decreasing redshift. The values found for Vrms 
are comparable to the level of kinematic disturbance measured in CiV systems (Ranch et 
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al. 2001a), suggesting that the metal deposition may have been recent in these systems, 
before the turbulence had time to dissipate. For a lower mean intergalactic metallicity of 
Z = 0.001 Zq, Vrms ~ 20 km s~^, which could result in significant broadening of the Lya 
absorption systems and help account for the discrepancy between the predicted and mea- 
sured line widths detected even in the optically thin systems (Meiksin, Bryan & Machacek 
2001). On the other hand, there would likely have been adequate time for the turbulence 
to dissipate by the time the winds could reach the underdense regions where the optically 
thin systems are expected to reside. Although the energy of the turbulence may be con- 
verted into heat as the turbulence dissipates, it is unclear how efficient this would be. Over 
a local Hubble time tniz), the energy associated with the turbulence, even if converted to 
heat, would also decay due to cosmological expansion. This is a topic that requires further 
exploration. 



5.2. Pressure Confinement of Bubbles 

All of the bubbles in our models are confined by external pressures if /* = 0.001, while 
all escape to the IGM if /* = 0.1. These results agree with our prediction in Figure 1 if 
the fraction of energy preserved in the bubbles v < 0.2. For detailed analysis of how ram 
pressure affects the evolution of the bubbles, we choose models with = 5 x 10^ Mq and 
/* = 0.01 at 2; = 3, 5, and 8, in which the efficiency of metal and mass ejection and of 
energy transport decrease steadily as a function of redshift. We also try to constrain u in 
equation (8) by comparing our results with Figure 1. 

We plot in Figure 12 the initial total pressure distributions as a function of radius in 
these three halos. We examine the sum of thermal and ram pressure because some kinetic 
energy is converted to thermal energy in shocks because of the artificial turn-off of cooling in 
the central regions, but the total energy is still conserved since we do not solve for coUisional 
coohng in the infall code (see §3.1). We also plot the total pressure distributions of bubbles 
at 45 degrees from the plane of the disks, measured at i = 85 Myr, for the halos at 2; = 3 
and 8. The ram pressure is included only for the non-colored gas moving with positive radial 
velocities. 

Figure 12 shows that the total pressure in the undisturbed halos in our models is a 
factor of five lower at r < i?^ than the ram pressure of the infalling halo gas calculated 
from equations (4) and (6) in §2, but the difference is less than a factor of two outside R^. 
Since pressure confinement mostly acts at r > i?„ for powerful bubbles with /* > 0.01, we 
can reasonably compare our results with the prediction made in §2. The external pressure 
increases as a function of redshift, consistent with the prediction from equation (4): pvf^ cx 
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(l + #. 

Now we try to constrain the bubble energy fraction u. All the bubbles in models with 
/* = 0.01 expand beyond Ryhy t — 50 Myr when the starbursts end, except along the planes 
of the galaxies where the ISM offers the most resistance (Fig. 5b). The bubbles generally 
grow to r = 2R-t, by t ^ 100 Myr, leaving most of the ISM and the fragmented shells of 
swept-up halo gas behind. The halo pressure has already begun to confine the bubble at 
z = 8 by t = 85 Myr (Fig. 12). Only a small fraction of original bubble energy remains 
at the swept-up shells in the form of kinetic energy. As a result, metal, mass, and energy 
feedback efficiencies are less than a few percent. Pressure confinement begins much later 
at 2; = 5, at t ~ 150 Myr, after the bubble grows beyond 27?^,. Since we assume that the 
low-density IGM begins at 2Ry, we compute moderately high ejection efficiencies of metals 
and halo mass of ~ 0.3-0.4 in this case. The pressure on the bubble at 2; = 3 is more 
than an order of magnitude less than at 2; = 8 (Fig. 12), so it expands beyond 2Ry without 
any significant infiuence from the infall, expelling most of the metals and returning a large 
fraction of infalling gas into the IGM. 

These halos fall just around the fine of pressure confinement with uf^ = 1 x 10~^ in 
Figure 1. We compute u in our models by summing the energy of interior gas with c < 0.6 
and swept-up gas with positive velocity, and dividing by Lmtf. More than half of the injected 
energy has already been lost in these halos by t = 50 Myr. Loss occurs by radiative cooling 
in the swept-up shells and in the mixed gas after the shells fragment. It continues as time 
progresses and is larger at higher redshift. Our resolution study shows that the dissipation 
of bubble energy seen in our simulations is physical, not dominated by numerical diffusion 
(see Fig. 11). 

We find the average energy fraction to be z/ ~ 0.2 between t = 50 and 100 Myr in the 
halos with Mh = 5 x 10^ Mq and = 0.01 at 2 = 3 and 5 and to he u < 0.2 at t = 50 Myr 
and z/ < 0.05 by t = 85 Myr in the same halo a.t z = 8. Therefore, we have u ~ 0.1-0.2 from 
our simulations for /* = 0.01, and so our results are consistent with the prediction made in 
Figure 1. We can also use equation (9) to define the confinement criterion in terms of virial 
temperature or circular velocity instead, finding that halos with > 1.0 x 10^ K or > 52 
km s~^ will begin to confine a starburst with = 0.01 and ~ 0.1. On the other hand, 
increases more than a factor of two to z/ > 0.6 in the same halos when /* = 0.1. This is 
because the cooling time of the swept-up shells Tsh becomes larger than the dynamical time 
of the bubbles within 2Ry, so the shocks remain adiabatic, and the fragmentation of the 
shells is suppressed. Bubbles even in galaxies with Vc ~ 280 km s~^ cannot be stopped by 
ram pressure if = 0.1 with u > 0.6. 

In summary, the ram pressure of the infalling halo gas strongly infiuences feedback 
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efficiency in dwarf galaxies, although it is primarily determined by the strength of star 
formation. The analytic results of §2 appear to offer a rehable guide to the importance of 
ram pressure. 

5.3. Fallback Timescales 

We compute the time for swept-up ISM and halo gas to cool and fall back to the center, 
which gives a minimum time before the next star burst. We define the accretion fraction as 

T{t) = Mi{R,t)/Mi{R,,0), (29) 

where Mi{R, t) is the gas mass with T < 10^ K within radius R at time t, and i?^, is the size 
of the galactic disk. Since we run our simulations for 200 Myr after the onset of starbursts 
or until the bubbles leave the grids, we first compute T(i) up to the end of the runs with 
R — i?*. 

If most of the gas is bound to the potential at that time, but hasn't yet fallen back to the 
center, we calculate T(t) for later times using a one-dimensional, ballistic approximation with 
R — Q. Bound gas parcels with v < v^sc are assumed to travel on radial, ballistic orbits in the 
potential, ignoring non-radial motions and acceleration of the accretion by cooling. We apply 
the approximation to gas within 2R^ at the time that feedback efficiencies are computed. 
This simple ballistic approximation was used by Zahnle & Mac Low (1995) to follow the 
ejecta of a typical Shoemaker-Levy 9 impact falling back onto Jupiter's atmosphere, and 
was found to give results consistent with observations. We also tested the validity of the 
method by computing the accretion fraction T(t) based on the ballistic approximation in a 
model in which we also compute T sirn{t) directly in the simulations. The bottom panel of 
Figure 13 shows that they agree within ~ 5%. We find similarly good agreement for other 
cases in which T simit) approaches unity so that we can make the comparison. 

Figure 13 also shows the accretion fractions T(t) as a function of time for three halos 
with M = 5 X 10^ Mq at ^ = 3, 5, and 8, with /, = 0.01 and 0.001. Because the ballistic 
model is only a one-dimensional approximation, it does not give exactly the same answer as 
the numerics at low values of T, as seen in Figure 13. We also define the fallback time t/ 
as the time when T(t/) = e~^, and give t/ for all our model halos in Table 3. The bump in 
T(t) for the z — % model occurs when ISM material falls back to the center almost at the 
free-fall velocity, creating accretion shocks. 

Since none of the mass escapes in models with /* = 0.001, T(t) quickly approaches unity 
within {).2tH{z). With /* = 0.01, feedback is more efficient in halos at lower redshift, but the 
timescale for gas to fall back remains < QAtH{z). This is rather significant in halos at 2; = 3 
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and 8 that have measured mass ejection efficiencies ^ ~ 1. Most disk gas appears tightly 
bound in galaxies with star formation efficiencies /* < 0.01. In such galaxies, subsequent 
starbursts will occur in pre-enriched gas with a fraction {1 — Cmet) of the metals from previous 
starbursts, yielding 



The presence of processed metals probably enhances star formation efficiency, which is the 
dominant factor influencing stellar feedback in our simulations. If the initial starbursts are 
very large, with /* = 0.1, the mass in the halos is blown away, and does not fall back within 
the local Hubble time tniz). Of course, real galaxies may never reach this star formation 
efficiency, as we discussed in §4. 

Our results suggest that it is difficult to keep most of the gas in a dwarf halo hot 
and diffuse and to prevent the central concentration of gas, once a dense disk-like object 
has formed. Star formation in small dwarf galaxies may occur in the form of repeating 
starbursts, although some other mechanism must explain the pauses longer than 10^ years 
between observed bursts of star formation. 

If feedback is to solve the overcooling and angular momentum problems in the formation 
of massive galaxies, it must prevent the initial collapse and concentration of gas into dwarf 
galaxies in the halos of the massive galaxies. Once much of the gas in the massive halo has 
cooled to form dwarf galaxies, subsequent starbursts in the dwarf galaxies have great difficulty 
heating and expelling it permanently. The dwarf galaxies then lose angular momentum from 
dynamical friction against the parent dark halo. Feedback from dwarf galaxies may prevent 
gas in nearby dwarf halos from collapsing, however, leaving a reservoir of diffuse gas to 
build a massive galaxy without the baryons losing too much angular momentum through 
dynamical friction. 



We study feedback from dwarf starburst galaxies with virial temperatures Ty > 10'^ K 
at high redshift, by numerically modeling the interaction of super bubbles driven by repeated 
supernova explosions with the ram pressure from cosmological infall. We compute ejection 
efficiencies for mass ^, metals ^met, and energy C, after the bubbles have expanded well 
beyond the virial radius, and estimate the timescale for disturbed gas to fall back and become 
available for subsequent starbursts. We have tried to consistently choose approximations that 
underestimate the effects of supernova feedback, so that our computations of the efficiency 
of mass, metal, and energy feedback and fallback time should be reasonably strong lower 
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6. Conclusion 
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limits to the actual values. We now list our conclusions. 

• Ejection and transport efficiencies are primarily determined by the efficiency of star 
formation With = 0.1, nearly all the metals produced in the starburst escapes 
iCmet ~ 1), niost of the disk ISM and large parts of the infalling halo arc blown away 

> 1), and accelerated halo gas carries the energy of the bursts out to the IGM 
(C > 0.5). With = 0.01, moderate to high feedback efficiencies are only observed in 
halos with Vc ~ 30-50 km s~^. With f^, = 0.001, none of the metals, the mass, nor the 
energy can escape the potentials. 

• The ram pressure of infalling halo gas can suppress the growth of starburst bubbles 
and thus prevent the metal, mass, and energy feedback, unless f^, > 0.1. The energy 
available to drive the bubble depends on both f^, and on the fraction of energy in 
the hot gas actually available to drive the bubble z/. The minimum halo mass for ram 
pressure confinement with a given uf^ can be specified by a cut-ofT in virial temperature 
T* or circular velocity v* at any redshift (Eq. 9). Our prediction is consistent with the 
values of u, ^met, ^iid C found in our simulations for a given With i/f^ — Ix 10"^, 
we can only expect high feedback efficiencies in halos with < T* 1 x 10^ K and 
Vc V* ^ 52 km s~^, above the limit for hydrogen line cooling, and below the Jeans 
mass in the reheated IGM. 

• We find that small dwarf galaxies with Vc ~ 30-50 km s~^ at 2; < 8 are efficient in 
enriching the IGM with metals, yielding Z > 10~^ Zq with = 0.01. The early metal 
enrichment may be able to explain the absence of kinematic disturbance observed in 
low-density Lyman a clouds. Larger dwarf galaxies with Vc ~ 100 km s^^ at z = 3-5 
might enrich the IGM more effectively, yielding Z ^ 10~^ Zq, because efficient cooling 
expected in such halos may enable /* ^ 0.01. The late metal enrichment is consistent 
with high metallicities and large turbulent motions observed in high-density Lyman 
a clouds. These galaxies form out of 1-2(7 peaks in our ACDM Gaussian density 
perturbation. 

• We expect mass ejection efficiencies ^ ~ 0.3-1 and energy transport efficiencies 0.1- 
0.3 from the dwarf galaxies with Vc ~ 30-50 km s~^ with = 0.01. However, most 
of the ISM remains bound unless /* = 0.1. Instead, a large amount of infalling gas 
is turned around by the galactic wind and escapes the potential. This outfiow may 
carry kinetic energy into the IGM, broadening Lya absorption systems, and to other 
nearby halos to prevent gravitational collapse of the gas in them. We estimate that 
the ensemble of outflows may provide turbulent support with typical Vrms ~ 20-40 
km s~^. Further study of this multiscale problem is required to determine whether 
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this mechanism of suppressing dwarf galaxy formation can solve the overcooling and 
angular momentum problems. 

• The timescale for swept-up mass to fall back to the center of a halo is short compared 
to the Hubble time, r/ < OAtniz) if /* < 0.01. We expect star formation efficiencies for 
subsequent starbursts to be larger than those of the initial starbursts that we studied, 
because ejected metals bound in the potentials will enhance gas cooling. Once gas in 
a halo cools to form a disk-hke object, it is difficult to blow away the disk gas. 

We must mention the caveats in our study: 

Our models assume an instantaneous starburst at the center of a galactic disk. This is 
obviously an oversimplification, and such instantaneous starbursts, especially with /* ^ 0.01 
are probably not realistic. We instead expect multiple OB associations or star clusters to 
form over a few hundred pc across the disk (e.g. Vacca 1996; Martin 1998). Each cluster 
will likely behave like our models however, so the net effect will be greater than we compute, 
consistent with our attempt to present lower limits to the effects of stellar feedback. Given 
the dominant influence of the star formation efficiency on the importance of stellar feedback, 
a self-consistent, physically motivated modeling of star formation will be required to further 
our understanding. 

We entirely neglected magnetic fields in this study. Magnetic fields can inhibit the 
formation of cold, dense shells, and suppress the fragmentation of the shells, perhaps reducing 
or preventing the dissipation of the bubble energy through the fragmented shells. Then 
the fraction of energy preserved in the bubbles u might be larger, giving higher feedback 
efficiencies in the dwarf galaxies that we modeled. On the other hand, magnetic pressure 
and tension might act to help confine the expanding bubbles (e.g. Tomisaka 1998), depending 
on the strength of the field. 

Our models are not three-dimensional. This has consequences for two main reasons. 
First, we compute the Rayleigh-Taylor instability in fragmenting shells during blowout in two 
dimensions. This has the effect of increasing the size of the fragments (Mac Low et al. 1989) 
because only ring-like modes of the instability can form. In three dimensions, the instability 
takes on a spike and bubble form, with the spikes representing the dense fragments. This 
does not change the gross topology, since even in two dimensions the shell still overturns and 
the wind can escape. However, a detailed study of the fragment properties cannot be done. 
As we restrict ourselves to integrated quantities, however, this probably does not strongly 
affect our results. Certainly the two models computed by Wada & Venkatesan (2003) in 
three dimensions appear to grossly agree with our results: a star formation efficiency of 
/* = 0.14 in a galaxy with 10^ Mq of gas leads to metal escape, while with = 0.014 the 
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metals remain confined. 

Second, our models assume spherically symmetric gas infall, and the presence of a 
well-formed, rotationally supported disk. It is also an oversimplification to compute metal, 
mass, and energy feedback efficiencies when the bubbles reach Rg — 2R„. Reahstically, 
a halo is embedded in the complicated weblike structures predicted by three-dimensional 
hydrodynamic simulations (e.g. Miralda-Escude et al. 1996; Zhang et al. 1998). Even inside 
its virial radius, gas accretion occurs along filaments, and the accretion shock remains highly 
aspherical even after the protogalaxy has collapsed to less than a few times the virial radius, 
(Abel et al. 1998; 2000; 2002). In addition, a rotationally supported disk may not form before 
a significant amount of star formation occurs. We note, however, that a more dynamical, 
filamentary background will create funnels through which superbubbles and galactic winds 
can freely expand, potentially increasing the efficiency of metal ejection and energy transport. 

Although these caveats identify important points requiring further investigation, none 
of them appear to invalidate our results; nor do they suggest that our lower limits are drastic 
underestimates of the actual effects of stellar feedback. 
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Table 1: Model galaxy properties 
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Table 2: Metal and mass ejection, and energy transport efficiencies 
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Table 3: Estimated fallback time r/ 
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Fig. 1. — We plot the minimum halo mass for ram-pressure confinement of bubbles as a 
function of redshift for 0.0001 < i//* < 0.005, where a fraction of injected energy preserved 
in the bubbles is z/ and star formation efficiency is The galaxies we choose to model 
(described below in section 4.1) are shown in filled diamonds. The starburst bubbles in all 
our model galaxies arc predicted to escape to intcrgalactic space if vf^ > 0.005 and to be 
confined by ram pressure of the infalling gas if uf^, < 0.0005. 
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Fig. 2. — We show the fraction of mass tracked by the tracer field with c < 0.6, Mcoior{t), over 
the total mass injected with c = 0, MsN{t), in our test simulations as a function of time up to 
10 Myr. The results are based on a standard model a.t z = 13, with halo background density 
Pbg = 200po, Lm — 10^° erg s~^, the primordial radiative cooling Aune and a resolution of 4 pc 
{thick solid line), the ClOO model with extra cooling, lOOAjj^e {solid line with stars), the CO 
model with no cooling {solid line with squares), the DIO model with the background density 
10 times higher than the standard model {dashed line), and the LR model with zones four 
times larger than the standard model {dotted line). The tracer field fails to pick up some of 
the injected gas due to the unresolved, high-density shells at the contact discontinuity. 
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z 

Fig. 3. — We choose seven second generation systems {filled diamonds) with virial temper- 
ature > 10^ K {thick solid line). As a comparison, the minimum halo mass in which the 
gas can cool through molecular hydrogen is shown in solid line, computed by Tegmark et al. 
(1997) and given by B. Ciardi. The Jeans mass after reionization is plotted in dash-dot-dot 
line. We also plot halo mass expected to collapse from ACDM density perturbations as la, 
2a, and 3a peaks {dashed line). 
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Fig. 4. — Density distributions for our model with halo mass Mh = 5 x 10^ M0 ai z = 8 
are shown at times 30 Myr (left), 50 Myr (middle), and 100 Myr (right), after the onsets of 
starbursts with star formation efficiencies /* = 0.1% (bottom), 1% (middle) and 10% (top). 
Corresponding mechanical luminosites are Lm = 1-8 x 10^^, 1.8 x 10^°, and 1.8 x 10^^ erg s~^ 
respectively. 
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Fig. 5. — Density distributions at t = 100 Myr for = 0.1% (left panels) and 1% (middle 
panels) and at t = 50 Myr for = 10% (right panel) for models with halo mass (a) 
Mh = 5 X 10^ M© at ;z = 13 and 8, (b) Mh = 5 x 10^ Mq at ^ = 8, 5, and 3, and (c) 
Mh = 5x 10^° M© at ;z = 5 and 3. 
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Fig. 6. — Density and temperature distributions of the bubble in the disk with M — 5 x 
10^ Mq at redshift S, &t t — 10 Myr after the onset of starburst with = 0.01. We plot 
them through the vertical z-direction from the galactic center. 
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r (kpc) 

Fig. 7. — Density distribution of model with halo mass — 5 x 10^ Mq at 2; = 8 right 
after the bubble blows out of the disk at i = 19 Myr with /* = 1%. The simulation was run 
with a resolution of 2.6 pc within the central 1 kpc^, decreasing smoothly to 26 pc beyond 
r z ^ 2 kpc. The cooled, swept-up shell of ambient gas fragments due to Rayleigh- Taylor 
instability, and mixes with the hot, low-density interior gas very effectively. 
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Fig. 8. — Figure shows the colored mass Mc divided by the total colored mass injected with 

c = 1 MsN as a function of e = 1 — c, for the halo with — 5 x 10^ Mq at 2; = 8 for 
star formation efficiencies = 0.001, 0.01, and 0.1. The masses M^. were measured when 
the bubbles reached 2Ry at t = 100 Myr when = 0.001 and 0.01 and at t = 50 Myr when 
= 0.1. The values of Mc are insensitive to the choices of e as long as they are very small, 
except for the case with = 0.1 because it was measured at earlier time. 
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Fig. 9. — The average metallicity of the universe Z from hales with mass between M^i^i^) 
and Mjnax at redshifts, 3, 5, 8, and 13 for ^metf* = 0.005. 
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t (Myr) 

Fig. 10. — Fractions of total injected energy escaping from the potential with Mh — 5 x 10^ 

Mq at ^ = 3: total energy escaping (dash- dot- dot line), total energy excluding thermal 
energy carried by swept-up gas liable to be radiated, defined as u (thick dash-dot-dot line); 
kinetic energy (solid line with stars) and thermal energy (dashed line with diamonds) carried 
by swept-up halo gas; and kinetic energy (solid line) and thermal energy (dashed line) carried 
by hot interior gas with tracer c < 0.6. 
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Fig. 11. — Resolution study of fractions of total injected energy remaining on the grid with 

M — 5x 10^ Mq at z = 8, with = 0.01. The main panel shows total energy (solid line with 
diamonds) and total energy excluding the thermal energy carried by swcpt-up gas (dashed 
line with stars) with standard resolution (5.2pc: thin lines) and doubled resolution (2.6pc: 
thick lines). The ratio of the results at the two resolutions is shown in the top panel. Note 
that the bubble with doubled resolution escapes the grid at t = 60 Myr. 
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Fig. 12. — Initial pressure profiles (sum of thermal and ram pressure) for halos with Mh = 
5 X 10^ Mq at z = 8 (thick solid line), 5 (thick dashed line), and 3 (thick dash-dot-dot line) 
as a function of radius r normalized by virial radius R^. We also plot pressure distributions 
of bubbles at 45 degrees from the plane of the disks in the halos at 2; = 8 (solid line) and 3 
(dash-dot-dot line), at t = 85 Myr after the onset of starbursts with /* = 0.01. The bubble 
at 2; = 8 has begun to be confined by external pressure, while the bubble at 2; = 3 freely 
expands to the IGM. 
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Fig. 13. — Accretion fraction T(t) for halos with Mh = 5 x 10^ Mq at z = 3 (dash- dot- dot 
line), 5 (dash line), and 8 (solid line) with star formation efficiency = 0.001 (above) and 
0.01 (below). The accretion fractions computed from the simulations are plotted in thick 
lines and those computed based on the ballistic approximation are plotted in thin lines. In 
the bottom plot, the ballistic approximation is also compared with the simulation directly 
for the z — 8 model, showing no more than 5% deviation. 



